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Research  into  methods  of  identifying  edges  in  a  noisy  scene 
has  been  an  active  field  of  investigation  for  many  years.  Treatment 
of  the  subject  may  be  found  in  many  books  written  over  the  past  de¬ 
cade  ([1],  [2],  and  [3])  and  many  different  approaches  are  proposed. 
Recently  a  survey  and  comparative  analysis  of  the  methods  was  made 

[7]. 

The  body  of  this  paper  is  segmented  into  4  parts.  In  the  first  we 
derive  and  define  a  "Moment  Operator"  which  we  show  to  work  well  for 
step  and  ramp  edges.  Then,  we  define  and  characterize  second  order 
edges  using  the  concept  of  the  rotation  of  a  point  in  a  vector  field 
and  develop  the  detector  analytically.  In  Section  3  we  develop  the 
algorithms  for  implementing  the  previously  defined  operators.  Final¬ 
ly,  in  Section  4,  these  algorithms  are  evaluated  using  ROC  curves  and 
compared  with  previously  known  techniques. 

The  detection  of  edges  to  isolate  objects  in  a  scene  is  motivated 
by  many  distinct  problems.  One  such  problem  arises  in  a  tracking  sys¬ 
tem  where  the  input  video  image  is  analyzed  and  the  object  to  be 
tracked  identified.  Subsequent  input  and  feedback  to  the  drive  con¬ 
trols  causes  the  sensor  to  re-orient  to  a  new  position  in  an  attempt 
to  maintain  the  same  x-y  coordinate  position  for  the  object  in  the 
field  of  view.  While  this  problem  motivated  the  research  that  led  to 
this  paper,  the  results  herein  discussed  are  much  broader  in  scope  and 
application.  The  constraints  imposed  by  this  problem  led  to  a  method 
that  is  useful  in  high  data  throughput  systems. 
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SECTION  1 .  EDGES  FROM  MOMENTS 


First  order  edge  detection  methods  work  in  the  following  way:  A  pic¬ 
ture  function  f(x,y)  is  transformed  to  another  picture  function  F(x,y) 
=  Tf(x,y)  in  such  a  way  that  the  edges  of  objects  in  the  scene  will  be 
in  the  set  ( (x,y) :F(x,y)>W}  for  some  W.  The  usual  method  is  to  trans¬ 
form  the  picture  using  T  equal  to  the  gradient  operator.  Different 
edge  detection  methods  correspond  to  different  numerical  approxima¬ 
tions  to  the  gradient. 

The  method  used  in  our  edge  detection  program  is  not  based  on  deriva¬ 
tives.  To  reduce  the  effect  of  noise,  this  edge  detection  method  uses 
integrals. 


The  reasoning  for  the  use  of  moments  to  find  edges  is  as  follows.  A 
digitized  picture  can  be  thought  as  a  lamina  whose  density  at  each 
point  is  f(x,y),  so  points  of  high  intensity  correspond  to  points  of 
high  density.  A  point  (a,b)  on  an  edge  in  the  original  function  (see 
Figure  1 )  would  correspond  to  a  point  in  this  lamina  (digitized  pic¬ 
ture)  with  high  densities  on  one  side  and  lower  densities  on  the  other 
side.  Thus  if  we  look  at  a  small  lamina  centered  at  point  (a,b)  and 
compute  the  center  of  mass  of  this  small  lamina,  we  can  expect  the 
center  of  mass  to  lie  within  an  area  of  high  densities. 


(1) 


(2) 


=  regions  of  high  density 


Figure  1.  Example  Center  of  Mass  Vectors  for  (1)  and  Edge  and 
(2)  a  Region  of  Uniform  Intensity. 


Suppose  we  now  look  at  a  point  (c,d)  such  that  the  densities  around  it 
are  fairly  constant.  Then  the  center  of  mass  of  a  small  lamina  about 
it  would  be  close  to  (c,d).  In  this  case,  a  vector  from  (c,d)  to  the 
center  of  mass  would  be  very  small  compared  to  a  vector  from  (a,b)  to 
the  center  of  mass  in  the  previous  case. 
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We  conclude  that  one  way  to  transform  f(x,y)  to  F(x,y)  such  that 
edges  of  the  original  picture  lie  in  the  set  F(x,y)>W  is  to  replace 
every  f(x,y)  by  the  length  of  the  vector  from  (x,y)  to  the  center  of 
mass  of  a  small  lamina  centered  about  (x,y).  That  is,  F(x,y)  is  the 
magnitude  of  the  vector  from  (x,y)  to  the  center  of  gravity  of  a 
square  lamina  centered  at  (x,y)  whose  density  is  given  by  the  picture 
function  f(x,y). 

Figure  2(b)  is  an  example  of  how  this  method  works  on  a  scene  (Figure 
2(a))  typical  of  those  we  study  at  WSMR. 

Once  the  coordinated  (X,Y)  of  the  center  of  mass  of  a  lamina  about 
(x,y)  are  calculated,  the  direction  of  the  edge  (if  any)  can  easily  be 
found .  Since  (X,Y)  points  to  where  the  intensity  of  the  picture  is 
the  highest,  the  direction  of  the__edge  is  perpendicular  to  the  direc¬ 
tion  of  the  vector  from  (x,y)  to  X,Y.  Ifj*e  take  (x,y)  =  (0,0),  then 
the  direction  of  the  edge  is  0  =  Arc  tan  (Y/X)  +  tt/2. 

Thus  this  model  gives  for  each  point  in  the  scene  a  quantity  that 
measures  the  probability  that  a  point  is  an  edge  point  and  a  direction 
which  is  the  direction  of  a  possible  edge  through  that  point. 

The  model  introduced  in  Section  I  will  not  work  for  roof  edges.  This 
is  because  at  the  ver^  peak  of  the  roof,  exactly  where  the  edge  is 
situated,  both  X  and  Y  are  equal  to  zero.  In  order  to  detect  roof 
edges  we  need  to  take  advantage  of  the  direction  information,  and  as 
Figures  6(a),  (b)  and  (c)  show  we  need  to  detect  the  shearing  cause 
by  the  change  in  direction  of  the  vector  field  at  the  edge  points. 

One  way  of  doing  this  is  by  using  a  tool  from  the  theory  of  vector 
fields,  namely  the  rotation  of  a  vector  field  about  a  point. 

SECTION  2.  SECOND  ORDER  EDGES 


After  a  scene  is  processed  by  the  moment  edge  detector,  each  point  is 
assigned  a  direction  and  a  magnitude.  In  effect  this  specifies  a  vec¬ 
tor  at  each  point  of  the  plan  in  question;  i.e.,  these  vectors  define 
a  vector  field  over  the  scene.  An  important  tool  in  the  study  of  vec¬ 
tor  fields  is  the  rotation  of  a  vector  field  (see  [4]  and  [5]).  To 
define  the  rotation  of  a  vector  field,  suppose  a  vector  of  the  vector 
field  4>  at  the  point  (x,y)  is  given  by 

*(x,y)  *  {<J>(x,y),  \p(x,y) } 

4>(x,y)  -  X(x,y) 

<Kx,y)  -  Y (x,y) 
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If  a  curve  T  on  the  plane  (scene)  is  given  in  the  form 

T:  x  =  x(t),  y  =  y(t)  a  <  t  <  b 

then  $(t)  *  {$[x(t),  y(t)],  ^[x(t),  y(t)]  is  defined  on  the  interval 
[a,b]  (see  Figure  3). 


Figure  3.  A  Curve  T  and  It's  Corresponding  Vector  Field  <I>(t) 

For  each  t  e[a,b]  there  is  determined  an  angle,  the  angle  in  radians 
between  $(t)  and  $>(a)  measured  from  $>(a)  to  $(t).  This  angle  is  a 
many  valued  function  of  t.  The  continuous  branch  of  this  function 
(vanishing  for  t  =  a)  is  designated  by  0(t)  and  called  an  angular 
function  of  the  field  $  on  a  curve  T.  The  rotation  of  the  field  $>  on 
the  curve  Y  is  defined  to  be 

[0(b)  -  0(a)] 


If  F  is  a  closed  Jordan  curve,  then  the  rotation  is  found  by  subdivid¬ 
ing  T  into  two  curves  (not  closed),  computing  the  rotation  of  each, 
and  adding.  In  the  following,  T  is  taken  to  be  a  small  circle  about  a 
point. 

We  can  write  the  rotation  as 
[0(b)  - 

With  0(t)  =  Arctan  Y/X  +  ir/2,  we  make  the  following  observations: 

(1)  If  0(t)  *  constant,  then  ■  0  and  y  *  0.  So  y  «  0  when 

x  *  a  point  on  the  edge  of  an  object  in  a  scene  (see  Figure  4)* 
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Figure  4.  Vector  Field  at  a  Step  or  Ramp  Edge  Point 

(2)  If  0  is  symmetric  about  x  and  T  is  a  small  circle  about  x  = 
edge  point  on  a  roof  edge,  see  Figure  5,  then  with  V  =  (where 

T1  -  one  half  of  the  circle  and  *  the  other  half) 

^ft^  C^t  =  +  *M0(t)  =  TT  +  IT  =  2  IT 

r  n  r2 


Figure  5.  Vector  Field  at  a  Roof  Edge  Point. 

An  example  of  how  these  observations  can  be  used  to  detect  second 
order  edges  appears  as  Figure  6(c)  and  (d). 
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SECTION  3.  ALGORITHMS  FOR  IMPLEMENTATION 


A.  Calculation  of  Moment 

Since  we  are  interested  in  real  time  applications  of  these  methods  we 
simplify  the  calculation  of  X  and  Y  by  setting 

M  *  /^h  /^k  f  (x  +  T,  y  +  u)  dt  du  =  1 

This  can  be  justified  by  observing  that  M/4hk  is  the  average  of  the 
intensities  over  a  small  neighborhood  of  (x,y)  and  so  this  value  can 
be  approximated  by  the  average  value  of  intensities  over  the  entire 
picture.  This  would  then  be  just  a  scale  factor  and  so  could  be  left 
out. 

To  calculate  the  integrals  involved  we  use  an  integral  formula  [6]  of 
order  0(h6).  The  formula  for  integration  is 
9 

/  /  F(x,y)  =  I  Wi*Di  with  W.  =  25/324,  W  =  10/81 
i  =  i  1 

and  if  we  apply  this  to  the  integrals  for  X  and  Y  and  factor  out  all 
scale  factors  we  get 

Y  =  5  *  (Di  -  D5)  +  4  *  (D8  +  D2  -  D6  -  Z>4) 

X  *  5  *  (D7  -  D3)  +  4  *  (D8  +  D6  -  D2  -  D4) 

and  use  abs  (X)2  -  abs  (Y)2  for  the  associated  magnitude.  If  we  sweep 
a  3  x  3  window  across  digitized  scene  D7  can  be  taken  as  the  upper 
left  hand  corner  while  D3  is  the  lower  right  hand  corner.  In  this 
case  the  direction  of  a  possible  edge  is  equal  to 

0  =  Arc  tan  (j^  —  3  +  tt/2 
Y  -  X 

B.  Calculation  of  the  Rotation 


The  vector  field  of  a  roof  edge  will  look  like  the  vector  field  of 
Figure  5.  So  to  find  roof  boundary  points  we  have  to  find  points  for 
which  in  a  neighborhood  of  such  a  point  /  d0  *  2tt 

The  smallest  region,  in  the  discrete  case  over  which  we  can  take  an 
integral  is  a  2  x  2  window.  Thus  our  algorithm  sweeps  a  2  x  2  window 
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across  a  scene  and  computes  the  integral  f  d0  for  each  of  these  win- 

c 

dows.  If  it  turns  out  that  this  integral  is  equal  to  2ir  then  those  4 
points  which  make  up  the  window  are  classified  as  boundary  points. 

To  calculate  the  integral  of  a  2  x  2  window  we  use  an  approximation 
4 

/de  =  Z  A  0 

X/ 

£  =  1 

computed  by  a  computer  program  given  in  Reference  [10]. 

For  the  purposes  of  this  experiment  the  procedure  used  to  generate  a 
file  which  is  the  file  of  detected  second  order  edges  in  the  follow¬ 
ing: 


1.  From  the  original  file  (scene)  two  files  are  generated;  one  (ACI) 
contains  SQRT  (X)2  +  (Y)2;  and  the  other  (ANG)  the  angle  of 

(0,  0  <_  0  <_  2555)  a  possible  edge. 

2.  From  the  ANG  and  ACI  files  one  new  file  AAA  is  created.  AAA  is 
created  by  sweeping  a  2  x  2  window  across  the  ANG  file.  The  rotation 
is  calculated,  and  if  a  point  is  classified  as  boundary  then  to  the 
corresponding  point  of  AAA  (initialized  at  zero)  is  added  the  average 
of  those  elements  of  ACI  that  have  the  same  subscripts  as  those  of  the 
2x2  window  being  swept  across  ANG. 

Examples  of  how  this  method  works  are  Figures  6(c)  and  2(d). 

SECTION  4.  EVALUATION 


The  methods  described  above  were  tested  on  disks  whose  edges  were 
step,  ramp  and  roof  edges.  The  step  and  ramp  edges  had  edge  height 
equal  to  16  while  the  roof  edge  was  constructed  by  beginning  at  the 
center  with  gray  value  equal  to  100  incrementing  by  one  to  gray  value 
equal  132  and  then  decrementing  by  one  to  gray  value  100.  All  files 
were  128  x  128  x  8. 

To  test  the  effectiveness  of  the  different  operations  considered  here 
we  added  Gaussian  noise  of  different  standard  deviation  to  achieve  a 
given  signal  to  noise  ratio  and  then  tested  the  algorithms  (Figure  7). 

16  2 

The  SNR  ratio  was  measured  in  db;  that  is,  we  used  SNR  =  10  log^^(— ) 

n 

where  a  =  standard  deviation  of  the  noise.  For  the  ramp  and  step 
m 

edges  we  used  SNR  -  4,  5,  6,  14  while  for  the  roof  edge  the 
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signal  to  noise  ratios  used  were  10,  11,  12,  . 20.  To  measure 
the  effectiveness  of  the  different  algorithms  we  graphed  PF  =  the 
probability  of  false  alarms  vs  PD  =  the  probability  of  detection 
(Figure  9(a)  and  9(b),  for  details  see  [71).  Figures  8(b)  and  8(c) 
contain  examples  of  processed  roof  edge  disks  with  SNR  =  13.  The 
graphs  of  PF  vs  PD  (ROC  Curves)  for  the  corresponding  operators 
appears  in  Figure  8(a). 

The  results  for  different  operators  and  step,  ramp  and  roof  edges 
appear  respectively  in  Figures  10(a),  10(b),  and  10(c).  These  graphs 
show  that  the  performance  of  the  moment  operator  is  in  all  cases  bet¬ 
ter  than  that  of  the  Sobel  operator.  A  significant  improvement  is 
obtained  by  first  applying  the  average  and  then  the  moment  operator. 
When  the  signal  to  noise  ratio  is  high  the  median  gives  better  results 
that  the  average,  but  there  is  a  crossover  point  at  which  the  average 
filter  gives  better  results  than  the  median. 
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